Hidden phonon highways promote photoinduced interlayer energy transfer in twisted transition metal dichalcogenide heterostructures

Vertically stacked van der Waals (vdW) heterostructures exhibit unique electronic, optical, and thermal properties that can be manipulated by twist-angle engineering. However, the weak phononic coupling at a bilayer interface imposes a fundamental thermal bottleneck for future two-dimensional devices. Using ultrafast electron diffraction, we directly investigated photoinduced nonequilibrium phonon dynamics in MoS2/WS2 at 4° twist angle and WSe2/MoSe2 heterobilayers with twist angles of 7°, 16°, and 25°. We identified an interlayer heat transfer channel with a characteristic timescale of ~20 picoseconds, about one order of magnitude faster than molecular dynamics simulations assuming initial intralayer thermalization. Atomistic calculations involving phonon-phonon scattering suggest that this process originates from the nonthermal phonon population following the initial interlayer charge transfer and scattering. Our findings present an avenue for thermal management in vdW heterostructures by tailoring nonequilibrium phonon populations.


Contents
1. Fitting models for time dependent MSD 2. Thermalized bilayer Molecular Dynamics simulations 3. First-principles phonon calculation using a perturbation theory approach 4. Supplementary figures and tables

Fitting Models for time dependent MSD
Two models can be used to represent the time-dependent response of the MSD.The first is a single exponential rise and double exponential decay model convolved with the instrument response function, where  is the amplitude of the response,  0 is the time of excitation,  is the instrumental resolution,  1 is the time constant for the fast exponential rise, and  3 is the time constant for the fast exponential decay due to vertical heat dissipation. 4 is the time constant for the slow exponential decay due to lateral heat dissipation, and  is the contribution of  4 to the decay in comparison with  3 . 4 is fixed to 6 ms as previously measured on WSe2/MoSe2 twisted bilayers on Si3N4 (48).Although our experimental time scale is not sensitive to such long time, it will add a non-zero offset to the end of the scan.This model well represents all monolayer MSDs and the MoS2 and MoSe2 layers of all heterobilayers.
For the WSe2 and WS2 layer of the heterobilayer, the model contained a fast exponential rise time constant ( 1 ), a slow exponential rise time constant ( 2 ), and the same fast exponential decay time constant ( 3 ) and slow decay time constant ( 4 ): where  is the contribution of  2 to the rise over  1 .
Since the time constants  are longer than experimental resolution , the model can be simplified and gives similar best-fit parameters: (3)

Mean square displacement and temperature
To quantitatively determine the Debye-Waller term, we perform molecular dynamics calculations to relate the temperature of individual MoSe2 and WSe2 monolayers to their in-plane mean-square displacements (MSD).To achieve this, we initialize both layers at the same temperature,  0  =  0  , from 40 K to 100 K, and subsequently calculate their final in-plane MSD at long timescales (1 ns).As shown in Fig. S4, the MSD displays the expected linear dependence on the temperature.Comparing Fig. S4 to the MSD peaks in Fig 2 .near 0.01 Å 2 suggests peak in-plane temperatures near 100 K for both the MoSe2 and WSe2 monolayers after the MoSe2 is optically pumped.We perform all MD simulations with the LAMMPS package incorporating both a Stillinger-Weber potential (63) for modeling interactions within a single TMD layer and Kolmogorov-Crespi (64) for capturing interlayer interactions, both of which are critical to capturing proper atomic relaxation in TMD bilayers.

Phonon eigenvector and dispersion
We use the PHONOPY package to compute, within the harmonic approximation, the phonon frequencies   and eigenvectors   () with crystal momentum wave vector  and branch index .These eigenvectors are obtained by diagonalizing the dynamical matrix constructed from the second-order force constants Φ  (,  ′ ′).The second-order force constants are defined as the change in force on atom ′ in unit cell ′ along the Cartesian direction  given a displacement in atom  in unit cell  along the Cartesian direction .We use a series of single-point calculations using the LAMMPS code to obtain the forces on all atoms in the system given some set of displacement of the lattice from equilibrium.This set of displacements is generated by PHONOPY to include the 0.01 Å displacements of each atom in all three Cartesian directions.The phonon dispersion is shown explicitly in Fig. S7.Phonon eigenvectors are projected onto layers such that each band is colored by the square of this projection.

Phonon lifetimes from phonon-phonon perturbation theory calculation
We use the PHONO3PY (55,65) package to calculate the phonon lifetimes of phonon modes in MoSe2/WSe2 bilayer heterostructures at several high-symmetry stacking positions and interfacial twist angles.The phonon lifetimes are determined from the imaginary part of the phonon selfenergy obtained from third-order force constants Φ  �,  ′  ′ ,  ′′  ′′ � which are defined by the change in force on atom  ′′ in unit cell ′′ along the Cartesian direction  given a displacement in atom  in unit cell  along the Cartesian direction  and a displacement in atom  ′ in unit cell  ′ along the Cartesian direction .These anharmonic, third-order terms are required to properly capture the three-phonon scattering process.To calculate the third-order force constant, we use PHONO3PY to generate an irreducible set of displacements of two atoms (,  ′  ′ ) with magnitude 0.03Å.We then perform many single-point calculations with LAMMPS to determine the change in force on the atom ( ′′  ′′ ) due to these displacements.However, this step quickly becomes computationally expensive for large lattices as the size of the displacement set scales as (3 atoms ) 2 and the size of the scattering matrix scales as (3 atoms ) 3 .To reduce the computational cost, we set a cutoff distance for the force.This sets the force between atoms beyond the cutoff to zero, reducing the number of force calculations to perform and the size of the scattering matrix in memory.We choose a cutoff distance of 7.5 Å, which accurately captures the intralayer interaction to second nearest-neighbor and the interlayer interaction to first nearest-neighbor.This approach accurately reproduces the phonon dispersion and phonon lifetimes of these bilayer TMD heterostructures.The generation of third-order force constants with force cutoff distances is performed by PHONO3PY.
Once constructed, we use the scattering matrix to compute the imaginary part of the self-energy which takes a form analogous to Fermi's golden rule, where Γ , () corresponds to the phonon linewidth of the phonon mode (, ),  is a phonon frequency,  is the Bose-Einstein occupation factor, and the phonon lifetime is calculated as 2Γ , � , � . (5) The construction and solution of the imaginary part of the self-energy from the third-order force constants are handled by the PHOEBE package (66).We calculate phonon lifetimes on a 10x10x1 grid for a 4x4 supercell of MoSe2/WSe2 bilayer heterostructures in the 0˚,  ℎ  stacking with 6 atoms in the primitive unit cell, and at finite interfacial twist angles 14˚ and 22˚ with 114 atoms and 42 atoms respectively.The lifetime for different phonons, resolved by their frequencies, is plotted in Fig. S8.Those phonons which are most likely to be produced in a non-thermal quantity during the initial hole relaxation after being generated by the pump are found in the shaded region.

Interfacial phonon-phonon scattering matrix elements in a layer-separable basis
The primary limitation of the above approach to quantify the non-thermal interfacial heat transfer is that the layer characteristics of the phonon modes involved in phonon-phonon scattering events are not well-defined.Naively, we might define the layer projection of a phonon in a bilayer system by the normalized displacement of the phonon modes within each layer: a phonon in layer, say, A would be a phonon whose normalized displacement is primarily in layer A. This would be an appropriate approach to calculate phonon-phonon interactions between phonon modes whose displacement is large, say, >99% in a single layer but fails when considering phonons that are nearly degenerate but very weakly hybridized due to the inherit degeneracy of the solutions.Any calculation of the interlayer heat transport involving an initial or final state that is layer-hybridized is therefore inaccurate and will underestimate the interlayer scattering time -and also will not be gauge independent for degenerate states.
To overcome this problem of characterization, we instead employ a layer-separable basis approach in which we rotate the phonon-phonon scattering matrix into a basis in which we can directly characterize all phonon modes by layer.We used a similar approach to Ouyang et al. (56), which calculated the phonon-phonon scattering elements by the off-diagonal components of the bilayer dynamical matrix rotated into a monolayer basis, but we extended this approach here to include the anharmonic effects from 3-phonon scattering processes.The separable monolayer basis, labeled below as AB, is constructed from the monolayer phonon modes, explicitly calculated on a single layer with the adjacent layer removed.These phonons are composed solely of displacements within a single layer and so are a good basis for isolating those phonon-phonon scattering events between layers.We then calculate the phonon-phonon scattering matrix elements as where F  (3) �,  ′ , "� is the third-order force constants in cartesian direction {, , } ∈ { �,  �, } calculated on the full bilayer system using LAMMPS classical force fields,  ,  are the eigenvectors of the interacting phonons in the monolayer basis with wave vector  and branch  with  = (, ) calculated on the separated monolayers using PHONOPY, and    are the monolayer phonon mode energies.In the separable monolayer basis, we can then easily extract those matrix elements which involves scattering from an initial phonon in layer A to two phonons, at least one of which is in the adjacent layer B. This selection procedure can be evaluated as Using the identical Fermi's golden rule approach as previously employed, we can then calculate the interfacial phonon lifetimes using these phonon-phonon scattering matrix elements with the interlayer phonon scattering selected.

Examining non-equilibrium phonon characteristics
In the main text, we discussed a representative fast interfacial phonon scattering process which originates from an acoustic  ≈  phonon in the MoSe2 layer.Our calculations show that the fast interlayer heat transfer time between WSe2 and MoSe2 originates from a  ≈  phonon with chiral displacement of the Se atoms about the xy-plane with a large out-of-plane displacement component throughout the oscillation period (Fig. S9).Previous time-domain reflectance experiments (TDTR) and DFT calculations show that the acoustic phonon modes with long-wavelength out-of-plane displacement carry a majority of heat between TMD layers (67).This particular phonon's large out-of-plane component could therefore similarly contribute to strong interlayer phonon coupling.

Lack of a substantial phonon highway in bilayers with different chalcogens
Our proposed mechanism depends on the presence of low-energy ZA phonons near the K point on both layers, and which are close in energy.Such a condition allows for the  ≈  phonon on the hot layer to be scattered to a  ≈  phonon on the cold layer with subsequent emission of a near-Γ, layer-hybridized acoustic phonon.The lowest-energy phonons at K primarily involve out-ofplane oscillation of the chalcogen atoms, and are hence nearly degenerate in bilayers with the same chalcogen atoms, such as bilayer WSe2/MoSe2, which we refer to as homochalcogen bilayer.
Similarly, we argue that this mechanism should be substantially reduced on a bilayer where the two chalcogen atoms are different, which we refer to as a heterochalcogen bilayer.
To test the hypothesis that the fast interlayer heat transfer should be weak on a heterochalcogen bilayer, we compute the phonon band structure of bilayer MoS2/WSe2 in Fig. S10.As expected, we find a much larger energy gap of 50 cm -1 between the two ZA acoustic bands at K at such a bilayer.We also compute the 3-phonon scattering processes in MoS2/WSe2 following the same approach as for MoSe2/WSe2 and find, as expected, longer scattering times associated with the interlayer phonon -about one order of magnitude larger in the heterochalcogen system compared to the homochalcogen bilayer (Fig. S11).We observe that there are still regions of fast interlayer scattering where the Mo acoustic bands overlap with the higher energy LA and optical bands of the W layer, consistent with our proposed mechanism.

MD simulation with nonthermal phonon distribution
For comparison to our perturbative theory approach, we perform an MD simulation on bilayer MoSe2/WSe2 with 0° twist angle to capture the temperature rise time of the WSe2 layer,   .However, instead of directly setting the temperature with a thermostat and thermalizing both layers, which generates a classical Maxwell-Boltzmann distribution of phonons, we add a non-thermal initial distribution of phonons to the MoSe2 layer by freezing in a number of specific phonons.If we initialize a molecular dynamics simulation with such a non-thermal distribution of phonons, the system will dynamically evolve and natively include all orders of phonon-phonon interactions, though assuming classical phonon statistics.
We introduce these phonons by displacing atoms from their position at time  = 0 according to the phonon eigenvectors.We add all the thermal energy in the simulation to ZA phonons close to K to simulate a highly non-equilibrium phonon population.The specific mode was chosen from those displaying fast interlayer scattering, as observed in Fig S6 .We evaluate the atomic displacement associated with a finite number of phonons as × 2Re� ⋅   , ()�, (8) where Δ  represents the displacement of an atom  having a mass of   and situated in the unit cell with a lattice vector of   .Here,   denotes the total number of unit cells within the supercell,  , () is the phonon polarization vector, normalized within the unit cell, with wave vector , band index , frequency   , and  , is the associated Bose-Einstein occupation at temperature T. In a supercell comprising 5400 atoms, we find that the effective number of acoustic phonons that need to be frozen in to increase the temperature of MoSe2 by 50K results in an average displacement of 0.05 Å of the Se atoms.This value was determined by post hoc analysis of the MD simulation to obtain the desired temperature using a linear scale factor mimicking a larger phonon occupation factor.
For the first mode, we choose a phonon from the ZA band near the K point, which displays a short interlayer scattering time.By freezing this phonon in the Mo layer at  = 0 and allowing the system to evolve under constant volume and energy, we observe the W layer temperature rise time of  MD W = 59 ± 5 ps (Fig. S10b).For completeness, we repeat this simulation for a phonon from the LA band at the K point, which shows a long interlayer scattering lifetime.In this case, we observe the W layer temperature rise time of  MD W = 246 ± 9 ps (Fig. S10c).These results should be compared to the thermalization time of  MD = 190 ± 7 ps reported in the main manuscript for a similar structure but for an initial thermal phonon within each layer (Fig. 3).

Supplementary Figures and Tables
Fig. S1.Example Sample Images.(a) 4˚ MoS2/WS2 and (b) 7˚ WSe2/MoSe2 on TEM grid.The heterobilayer on the Si3N4 membrane window is the grey square in the center.The surrounding area is the Si grid.For each sample, the heterobilayer homogenously covers the entire window and exhibits minimal cracks or multilayer regions.

Fig. S3 .
Fig. S3.Debye Waller Response for MoS 2 /WS 2 .(a)(b) Integrated Bragg peak intensities as a function of pump-probe delay time for Bragg peaks in the MoS2 and WS2 layers of a 4° MoS2/WS2 heterobilayer pumped on resonance with the MoS2 A exciton with a fluence of 1 mJ/cm 2 at 50K.(c) Log intensity change for 7 orders of Bragg peaks at delay time of 12 ps, plotted against the reciprocal lattice vector of the Bragg peak squared, with the corresponding linear fit of the Debye-Waller model.

Fig. S4 .
Fig. S4.MSD and Temperature.MSD of the MoSe2 and WSe2 layers in a 7° MoSe2/WSe2 heterobilayer as a function of temperature.Fitted line shows expected linear relationship between MSD and temperature T.

Fig. S6. Three-phonon scattering lifetimes for MoSe 2 /WSe 2 .
Fig. S6.Three-phonon scattering lifetimes for MoSe 2 /WSe 2 .Interlayer three-phonon scattering lifetime computed for an initial phonon of primarily MoSe2 character.The phonon dispersion of the interacting system (thick lines) is also compared with that of an isolated layer of MoSe2 (faint red dots) and WSe2 (faint blue dots).Note that phonon dispersion associated with an isolated layer of MoSe2 is nearly completely occluded by that of the interacting bilayer system.The fastest

Fig. S7. MoSe 2 /
Fig. S7.MoSe 2 /WSe 2 calculated phonon dispersion.Phonon dispersion for the MoSe2/WSe2 heterostructure projected onto the respective layers, computed with PHONOPY (62).Layer hybridization is minimum apart from band crossings and long wavelength acoustic modes.Note that the lowest-energy phonon bands near K are not strongly hybridized.

Fig. S8 .
Fig. S8.Phonon lifetimes at different twist-angles.Phonon lifetimes at 300K for 0˚, 14˚, and 22˚ twisted MoSe2/WSe2 heterobilayers.Phonon lifetimes are calculated from the scattering matrix triplet elements between an initial phonon mode and two scattered phonon modes which preserve energy and crystal momentum.Grey region highlights those initial phonon modes which originate from acoustic modes near the K point of the BZ.Acoustic phonons between energies of 3.5 THz and 5 THz are primarily produced during the K-Γ hole scattering of the intra-to-interlayer MoSe2/WSe2.Minimum and maximum phonon lifetimes for this region is 3.5 ps and 35.1 ps, respectively.

Fig. S9 .
Fig. S9.Dynamics of a MoSe 2 phonon with fast interlayer phonon-phonon scattering.A particular  ≈  wavevector phonon of MoSe2 along the K-M line which shows fast interlayer phonon-phonon scattering.The phonon has an energy of 17.4 meV (period of ~0.24 ps).The primary characteristic of this phonon is a chiral motion of the Se atoms along the in-plane direction and a displacement of the transition-metal atom along the out-of-plane direction.The phonon eigenvector was calculated using PHONOPY(62).

Fig. S11. Three-phonon scattering lifetimes for MoS 2 /WSe 2 .
Fig. S11.Three-phonon scattering lifetimes for MoS 2 /WSe 2 .Interlayer three-phonon scattering lifetime computed for an initial phonon of primarily MoS2 character in a MoS2/WSe2 heterochalcogen bilayer.The interlayer phonon lifetime of the initial Mo phonon (thick lines) is overlayed on the phonon dispersion of an isolated layer of MoSe2 (faint red dots).The phonon dispersion of an isolated layer of WSe2 (faint blue dots) is shown for comparison.The fastest interlayer scattering regions occur near the band crossings of MoS2 and WSe2 phonons, which are dramatically reduced along the ZA band compared to the MoSe2/WSe2 due to the presence of the energy gap between monolayer energies.The vertical red lines delimitate the expected range of initial phonon wavevectors originating from the hole relaxation in the MoS2 layer.

Fig. S12 :
Fig. S12: Interlayer thermalization time from MD simulations.a) Interlayer thermalization time from MD simulations for the homochalcogen bilayer MoSe2/WSe2 with a 0˚ twist angle on a commensurate R X H stacking.We chose two phonons to analyze: one from the ZA band with a short interlayer scattering time (green) and one from the LA band with a long interlayer scattering time (pink).b) Interfacial thermalization for a nonequilibrium phonon distribution.The MoSe2 and WSe2 layers are initially equilibrated with an NVT thermostat to 50 K and 50 K, respectively.At t = 0, a specific phonon mode from the ZA band close to K with a short interlayer scattering time is frozen into the Mo layer.The effective temperature of the layer is raised to 100K due to the nonequilibrium population.Immediately after, the temperature constraints are released, and the W layer temperature rises with a time scale τ MD = 59 ± 5 ps, considerably faster than that found in the thermally distributed simulation.c) Similar to (b), but we initially freeze the highlighted phonon mode from the LA band with long interlayer scattering time.When the temperature constraints are released, the W layer temperature rises with a time scale τ MD = 246 ± 9 ps.

Table S1 . Reciprocal lattice vectors used for Debye-Waller model.
Different Bragg orders and the corresponding Q 2 hk0 values for each monolayer.Lattice constants of 3.18 Å used for MoS2 and WS2 and 3.32 Å used for MoSe2 and WSe2 (68).